

#define RND_BUFFER	16777216
//#define UINT32_MAX (0xffffffff) //this is redefined on OSX?
#define SKEW_SEEDS	256
uint64_t u01n;
uint32_t *rnd_array=NULL;
sfmt_t sfmt;;


uint32_t *skew_seed(bool random) {
uint32_t *seed=new uint32_t[SKEW_SEEDS];
time_t tick;
uint j=0;
for (uint i=0;i<SKEW_SEEDS;i++) {
	if (random) {
        tick=clock();
        for(j=0;(uint)clock()==tick;j++)
                ;
	} 

        seed[i]= j%256;// only use lower 8 bits..
}
return seed;
}

float u01() {
	u01n++;

	if (u01n%RND_BUFFER==1) {
		if (rnd_array==NULL) {
			rnd_array=new uint32_t[RND_BUFFER];
			uint32_t *seeds=NULL;
			sfmt_init_by_array(&sfmt, seeds=skew_seed(RANDOM_SEED), SKEW_SEEDS);
			free(seeds);
		}
		sfmt_fill_array32(&sfmt, rnd_array, RND_BUFFER);
		u01n=1;
	}

return  rnd_array[u01n-1]/(float)UINT32_MAX;

}

float n01(float u1, float u2, bool use_sin) {

if (use_sin)
        return sqrtf(-2*logf(u1))*sin(2*PI*u2);
else
        return sqrtf(-2*logf(u1))*cos(2*PI*u2);

}

float simple_n01() {

return n01(u01(),u01(),0);

}

float chi2 (uint n) { //probably won't be needing this...

	float sum=0,u1=0,u2=0;
	for (uint i=0;i<n;i++) 
		if (i%2==0) {
			u1=u01(), u2=u01();
			sum+=powf(n01(u1,u2,0),2);
		} else
			sum+=powf(n01(u1,u2,1),2);
	
	return sum;
}


